Transition to an Insulating Phase Induced by Attractive Interactions 
in the Disordered Three-Dimensional Hubbard Model 
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We study numerically the interplay of disorder and attractive interactions for spin- 1/2 fermions 
in the three-dimensional Hubbard model. The results obtained by projector quantum Monte Carlo 
simulations show that at moderate disorder, increasing the attractive interaction leads to a transition 
from delocalized superconducting states to the insulating phase of localized pairs. This transition 
takes place well within the metallic phase of the single-particle Anderson model. 



PACS numbers: 71.30.+h, 74.20.-z, 71.10.-w 



Two limiting cases of the non-trivial problem of quan- 
tum transport in three dimensions (3D) in the presence 
of disorder and attractive interactions between particles 
were worked out by Anderson in the late 1950s In 
the limit of weak interactions, the increase of disorder 
leads to the Anderson transition from metallic transport 
to the localized insulating phase [jD. Whereas, in the 
absence of disorder, attractive interactions between spin- 
1/2 fermions lead to the appearance of BCS superconduc- 
tivity which is not affected by the introduction of weak 
disorder [Q. These limiting cases have been extensively 
investigated and detailed information is now available 
for the one-particle Anderson transition (see e.g. (| ||]) 
and the weakly disordered BCS superconductor (see e.g. 
^-||). However, a theoretical treatment of the intermedi- 
ate regime, where both disorder and interactions are im- 
portant, is difficult due to the absence of relevant small 
parameters. New results on the physical properties of 
transport in this regime are therefore of fundamental in- 
terest. Furthermore, an understanding of this realistic 
regime would contribute to the interpretation of recent 
experiments on 3D superconductors, where both disor- 
der and interactions are naturally present in the physical 
systems studied P-|l2|. Indeed, the experimental results 
of Ref. H indicate an intriguing correlation between the 
Anderson transition in a strong magnetic field and opti- 
mal doping for the superconducting transition tempera- 
ture. Also, an explanation of the unusual resistivity de- 
pendence on magnetic field observed in Ref. |i"l| ] requires 
a better understanding of the interplay between disorder 
and attractive interactions. In addition, recent break- 
throughs in cold-atom experimental techniques have pro- 
vided new possibilities for investigations of interacting 
atoms on 3D optical lattices, leading to the observation 
of a superfluid to Mott insulator quantum phase tran- 
sition for ultra-cold atoms [Q. These extremely clean 
experiments open unprecedented possibilities for precise 
studies of lattice models with experimentally tunable in- 
teractions and provide new challenges for theoretical in- 



vestigations. 

Numerical simulations provide a valuable tool for the 
study of the non-trivial regime where both interactions 
and disorder play a relevant role. Among various nu- 
merical approaches, quantum Monte Carlo methods con- 
stitute the most promising possibility for the simula- 
tion of 3D systems with a large number of particles 
|]i~4] [ijj . These methods have several advantages com- 
pared to exact-diagonalization approaches which are re- 
stricted to a relatively small number of particles H,ff9[ . 
In this work we use the projector quantum Monte Carlo 
method (PQMC) to perform numerical simulations of the 
3D Anderson transition in the presence of attractive in- 
teractions. To the best of our knowledge, this approach 
allows us to explore this new regime for the first time. 

To investigate the interplay of disorder and attrac- 
tive interactions, we study numerically the disordered 3D 
Hubbard model with N fermions on a cubic lattice with 
L 3 sites. The Hamiltonian is defined by 
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where c\ a (ci CT ) creates (destroys) a spin-1/2 fermion at 
site i = (i x ,i y ,i z ) with spin a, rii, 



is the cor- 
responding occupation number operator. The hopping 
term t between nearest neighbor lattice sites parameter- 
izes the kinetic energy and the random site energies ei are 
homogeneously distributed in the interval [—W/2,W/2], 
where W determines the disorder strength. The param- 
eter U measures the strength of the screened attractive 
Hubbard interaction (U < 0) and periodic boundary con- 
ditions are taken in all directions. At U — 0, the Hamilto- 
nian (Q) reduces to the one-body Anderson model, which 
exhibits a metal-insulator transition in three dimensions 
p-pl. For W = the Hamiltonian corresponds to the 
clean attractive Hubbard model, with a superconducting 
ground state in 3D. 

We study this model by the PQMC method, which is 
an efficient method for the investigation of the ground 
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state properties of interacting fermion systems fll5| , f20|] . 
For attractive Hubbard interactions (U < 0), there is no 
sign problem and the method is exact up to discrete time- 
step and statistical errors, which can be well-controlled. 
We consider N = 14, 32, 62, and 108 particles on a cubic 
lattice of size L = 3,4,5, and 6, respectively, at an ap- 
proximately constant filling factor v = N/(2L 2 ) « 1/4, 
2 < W/t < 10, and U/t = -4. For each disorder real- 
ization, we used a discrete Trotter decomposition with a 
time step At = 0.1 and projected through 60 time steps. 
In total we carried out 3000 Monte Carlo sweeps for each 
simulation, with approximately 1000 sweeps for equili- 
bration. With these PQMC parameters we obtained 
good convergence of the computed physical quantities. 
The results are averaged over Nr > 12 disorder realiza- 
tions, except for the most time-consuming simulations at 
L = 6, where Nr = 6. The simulations are carried out in 
the sector with the total spin component S z = 0. This is 
sufficient to study the ground state properties, since all 
of the eigenvalues of the Hamiltonian ([!]) belong to the 
spectrum of the S z = subspace. 

A quantitative measure of the localization properties 
of the system can be obtained, even in the presence of 
interactions. It is based on the probability density dis- 
tribution for an added pair at the Fermi edge. This dis- 
tribution is approximately equal to the charge density 
difference Sp(i) = p(i,N + 2) — p(i,N), where p(i) is 
the ground state charge density at site i (%2i5p(i) = 2). 
The values of p(i, N), p(i, N + 2) are obtained from two 
independent PQMC simulations for the same disorder re- 
alization. For U = 0, this difference is identical to the 
single-particle probability distribution for the eigenstate 
iPf{i) of the Anderson model at the Fermi level, with 
Sp(i)/2 — |V>F(i)| 2 - From this distribution, we obtain 
the inverse participation ratio (IPR) for an added pair, 
£ = (^;(#/?(i)/2) 2 ) _1 , which determines the number of 
sites visited by this pair. At U = 0, this quantity re- 
duces to the usual one-particle IPR at the Fermi level. 
For U ^ 0, since the interaction is screened and short- 
ranged, £ still determines the localization properties of 
pairs in the vicinity of the Fermi level. 

A typical example of Sp(i) for an added pair is shown in 
Fig.[j] for a single disorder realization. For graphical rep- 
resentation 5p(i) is projected on the (x,y) plane, giving 
5p p (i x ,iy) = J2i z Sp(i x ,i y ,i z ). The plots of Fig|] (left) 
refer to the weakly disordered regime (W/t = 2). They 
clearly show derealization of the added pair, both for 
free particles (top, U = 0) and for the attractive Hub- 
bard model (middle, U/t = —4). It should be noted that 
this disorder strength is much smaller than the critical 
disorder strength of the 3D Anderson transition which 
takes place at W < W C {U = 0) w 16.5* |l|. Therefore, 
these results confirm Anderson's theorem pj according to 
which the Cooper pairs remain delocalized at weak dis- 
order. Thus, the superconducting phase is not affected 
by weak disorder, since Cooper pairs can be formed by 



pairing the time-reversed eigenstates of the correspond- 
ing non-interacting disordered problem. 




FIG. 1. Distribution of charge density difference for an 
added pair, Sp p , projected on the (x,y) plane for a 6 x 6 x 6 
lattice for the same single disorder realization, with W/t — 2 
(left) and W/t = 7 (right), N = 108. Top: exact computation 
for U = 0, £ = 70; 55 (left; right). Middle: PQMC calculation 
for U/t = -4, £ = 48; 6.5 (left; right). Bottom: BdG mean 
field calculation for U/t = -4, £ = 132; 25 (left; right). 

A qualitatively new situation appears at a stronger dis- 
order strength W/t = 7 (Fig.|l| (right)). Here, the two 
added particles are delocalized at U = 0, since we are still 
inside the metallic single particle phase W < W C (U = 0). 
On the contrary, a pronounced peak appears for 5p(i) 
at U/t — —4 clearly showing the formation of localized 
pairs. This is borne out by the IPR which drops from 
£ = 55 at U = to £ = 6.5 at U/t = -4. This effect 
gives an indication that attractive interactions induce lo- 
calization in the metallic regime of the non-interacting 
model, leading to the formation of bi-particle localized 
states. 

The results of the Bogoliubov-de Gennes (BdG) mean 
field calculation J6|,^2| , are shown in the two bottom plots 
of Fig.|j]. We note that, within BdG approach, at weak 
disorder strengths (W/t = 2), interactions smooth out 
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charge fluctuations leading to an increase of the 1PR, 
from £ = 70 at U = to £ = 132 at U/t = -4. On 
the contrary, at stronger disorder strengths (W/t = 7) 
interactions slightly favor localization even within BdG 
approximation where the IPR drops from 55 to 25 when 
U/t goes from to —4. This happens because the mean 
field treatment of interactions introduces a site depen- 
dent Hartree shift U H {i) = \U\p(i)/2 [§§]. At strong dis- 
order, when the charge density p(i) is highly inhomoge- 
neous, this term acts as an additional disorder potential 
p3| . However, the main problem with the BdG approach 
is that many local minima appear and convergence be- 
comes problematic even at moderate disorder strengths 
and system sizes (W/t > 7,L > 6). Furthermore, it is 
clear from the comparison with the PQMC charge den- 
sity difference at W/t = 7 that important effects beyond 
mean field cannot be reproduced within BdG approach. 
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FIG. 2. Inverse participation ratio < £ > averaged over 
disorder realizations, as a function of disorder strength W for 
a 6 x 6 x 6 lattice, at U — (empty circles) and U/t = —4 
(filled circles). Dotted lines show linear fits to the data, the 
dashed line represents £ = 1 (see text) and error bars indicate 
statistical errors. 

A more quantitative description of the localization ef- 
fect induced by attractive interactions can be seen from 
the dependence of IPR < £ > on disorder strength, W, 
shown in Fig. 2, with < £ > averaged over realizations 
of disorder. The data clearly show that interactions lead 
to a significant reduction of < £ >. For {7 = 0, the 
dependence of < £ > on W is characterized by two dis- 
tinct regions: a relatively flat region for large W where 

< £ > slowly approaches the asymptotic value of 1 and 
another region in which < £ > grows with decreasing 
W. At very weak disorder < £ > remains bounded by 
the total number of lattice sites. The linear fit of the 

< £ > data in the second region crosses the < £ >= 1 
line at W C (U = 0) « 16t which is very close to the 
exact value of W for the Anderson transition for non- 



interacting particles, W c w 16. 5t 0). A similar analysis 
carried out for the < £ > data in the presence of interac- 
tions gives transition from delocalized to localized pairs 
at W C (U = -At) w 9t. 




FIG. 3. Dependence of the scaled inverse participation ra- 
tio < £ > /L 2 on the linear dimension of the system L, for 
U = (top) and U/t = -4 (bottom), with W/t = 3 (trian- 
gles), 5 (diamonds), 7 (squares), and 10 (circles). Error bars 
show statistical errors. The straight line fits show the aver- 
age dependence on L for the extremal values of W studied, 
designating the transition from superconducting to insulating 
behavior for U/t — —4. 

Further evidence for the interaction induced transition 
comes from a finite-size scaling analysis. The relevant 
dimensionless quantity for such an analysis is the sys- 
tem conductance g(L) pp^|. For large g(L), the macro- 
scopic transport theory in the 3D delocalized phase gives 
g(L) cx L cx £/L 2 since £ oc L 3 for delocalized wave func- 
tions. For the localized phase, £ is determined by the lo- 
calization length I, and is independent of the system size 
L, with £ ~ I 3 . Hence, the ratio £/L 2 falls off as 1/L 2 
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in this regime. Therefore the transition point g(L) w 1 
can be located from the condition £(L)/L 2 = const. The 
results of the finite-size scaling analysis for the scaled ra- 
tio £/ ' L 2 are shown in Fig.|| for 3 < L < 6. The range 
of disorder values studied (3t < W < lOt) corresponds 
to the metallic side of the single-particle Anderson tran- 
sition. Therefore the scaling analysis at U = shows 
an increase of £/L 2 with the system size. A strikingly 
different situation appears at U — —At: at W/t = 3 and 
5 the scaling ratio £,/L 2 still grows with the system size, 
while at W/t = 7 and 10 this quantity drops with L. 
This suggests the appearance of a superconductor to in- 
sulator transition induced by attractive interactions, with 
the transition point W C (U — —At) k, 6t at the thermody- 
namic limit. This value is in reasonable agreement with 
the value obtained for a single system size in Fig. 2. 
A precise location of the transition point would require 
a significant increase of the system sizes and a larger 
number of disorder realizations. However the results ob- 
tained in the present study clearly show the transition to 
an insulating phase at disorder strengths being less than 
half the value of the critical disorder strength for the 
single-particle Anderson transition. Hence, in the pres- 
ence of attractive interactions, the insulating phase pen- 
etrates inside the metallic non-interacting phase. This 
unexpected result can be understood on the basis of the 
following physical argument p5[ . The attractive interac- 
tion creates pairs of effective mass to* twice as large as 
the single fermion mass to. This halves the effective hop- 
ping term t* oc I/to* = 1/2to which induces an increased 
effective disorder W/t* — 2W/t and thus enhances local- 
ization effects. Such an argument predicts the decrease 
in critical disorder strength by a factor of 2, which is in 
reasonable agreement with our results. 

In conclusion we show that in disordered systems, at- 
tractive interactions that lead to superconductivity at 
weak disorder also lead to the insulating phase of local- 
ized pairs at moderate disorder strengths, well within the 
metallic phase of non-interacting fermions. Thus, by in- 
creasing the attraction strength, one can go from super- 
conducting to insulating behavior. This is also possible 
by increasing disorder. Of course, experimentally, it is 
not easy to vary the interaction and disorder strengths. 
However, indirectly, this can be achieved by introducing 
a relatively strong magnetic field. This magnetic field can 
increase the effective disorder strengths since it forces an 
electron to return to the impurity |p6| . At the same time, 
it also effectively decreases attraction by pair breaking. 
Thus, the increase of magnetic field may first increase 
the disorder and drive the system from superconductor 
to insulator with localized pairs as seen experimentally 
in fnj]. Further increase of magnetic field breaks pairs 
and leads to a transition from insulating phase of local- 
ized pairs to a metallic phase of almost non-interacting 
fermions. Such a scenario leads to a transition from su- 
perconductor to insulator, followed by an insulator to 



metal transition with increasing magnetic field, in quali- 
tative agreement with experimental observations jnj. 
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